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SUMMARY 


A  new  method  is  presented  for  calculating  the  scattering  of  an  arbitrary 
electromagnetic  wave  by  a  bounded,  perfectly  conducting  body  of  general  shape. 
The  strategy  is  to  replace  the  corresponding  exterior  boundary-value  problem  by 
an  approximate  problem  on  the  boundary  of  the  scattering  body.  This  Involves 
the  introduction  of  a  certain  bilinear  form  and  non-local  boundary  operator, 
together  with  the  use  of  a  special  class  of  known  solutions  of  the  reduced 
Maxwell's  equations  satisfying  the  Sommerfeld  radiation  boundary  conditions  at 
Infinity.  Two  computer  programs  implementing  this  method  are  described  and 
numerical  results  showing  the  successful  application  of  this  method  to  some 
model  problems  are  presented. 

B.  P.  DESATAGE 
By  direction 
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1.  INTRODUCTION 


The  numerical  calculation  of  scattering  by  an  arbitrarily  shaped  perfect 
conductor  or  dielectric  body  has  received  much  attention  in  recent  years.  The 
problem  is  important  in  many  areas  including  the  design  of  waveguides  and  missile 
fuzes,  the  assessment  of  damage  by  an  electromagnetic  pulse,  and  the  study  of  the 
biological  effects  of  microwave  radiation.  Mathematically,  the  problem  has  the 
form  of  an  exterior  boundary  value  problem  for  a  system  of  elliptic  partial 
differential  equations.  The  problem  has  particular  difficulties  arising  from  the 
large  number  of  unknowns,  the  unbounded  domain,  and  the  fact  that  the  solution  is 
oscillatory  for  large  values  of  frequency. 

The  principal  techniques  for  the  numerical  solution  of  scattering  problems 
are  the  method  of  Integral  equations,  and  the  finite  difference  or  finite  element 
method.  We  discuss  briefly  these  methods. 

In  the  method  of  integral  equations,  one  starts  with  a  formulation  of  the 
boundary  value  problem  as  an  integral  equation  on  the  surface  of  the  scattering 
body;  the  unknown  may,  for  example,  be  the  current  density  on  the  surface.  An 
excellent  survey  of  such  formulations,  in  the  case  of  the  reduced  wave  equation, 
is  given  in  [1].  The  integral  equation  is  solved  numerically,  either  by  using  a 
quadrature  formula  or  a  Galerkin  method  to  reduce  the  problem  to  a  finite  system 
of  linear  equations.  The  principal  advantage  of  the  method  is  that  the  dimension 
of  the  problem  has  been  reduced;  one  must  find  a  vector  field  on  the  two  dimensional 
surface  of  the  scattering  body,  Instead  of  in  the  three  dimensional  exterior 
region.  On  the  other  hand,  the  system  of  linear  equations  has  a  full  (non  sparse) 
coefficient  matrix.  Also,  the  integral  equation  contains  a  weakly  singular 
kernel  that  comes  from  the  fundamental  solution  of  the  problem;  the  resulting 
surface  Integrals  may  be  difficult  to  evaluate  with  sufficient  accuracy.  Examples 
of  the  integral  equation  approach  are  contained  in  [2,3]. 

With  the  second  method,  the  boundary  condition  at  infinity  is  replaced  by  a 
boundary  condition  on  the  surface  of  a  large  sphere  containing  the  scattering 
body.  This  approach  is  used  in  [4].  The  resulting  boundary  value  problem  is 
discretized  by  means  of  a  finite  difference  or  finite  element  method.  This 
method  has  the  advantage  of  producing  a  simple,  sparse  coefficient  matrix. 

On  the  other  hand,  the  order  of  the  matrix  will  be  larger,  because  the 

unknown  is  now  a  vector  field  in  a  three  dimensional  domain.  A  judicious  use  of 


1.  Kleinman,  R.  E.  and  Roach,  G.  F. ,  "Boundary  Integral  Equations  for  the  Three 
Dimensional  Helmholtz  Equation,"  SIAM  Rev. ,  Vol.  16,  pp.  214-236,  1974. 

2.  McDonald,  B.,  Friedman,  M. ,  and  Wexler,  A.,  "Variational  Solution  of  Integral 
Equations,  "  IEEE  Transactions  on  Microwave  Theory  and  Technique.  Vol.  MTT-24, 
pp.  237-248,  1974. 

3.  Livesay,  D-  E.  and  Chen,  K-M. ,  "Electromagnetic  Fields  Induced  Inside 
Arbitrarily  Shaped  Biological  Bodies,"  Ibid.,  pp.  1273-1280. 

4.  Baylis,  A.,  Gunzburger,  M. ,  and  Turkel,  E. ,  "Boundary  Conditions  for  the 
Numerical  Solution  of  Elliptic  Equations  in  Exterior  Regions,"  Report  80-1, 
ICASE,  NASA-Langley,  Hampton,  Va. ,  1980. 
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graded  meshes  ac  a  large  distance  from  the  scattering  body  will  help  alleviate 
this  problem  [5].  Examples  of  the  use  of  finite  differences  or  finite  elements 
are  contained  in  [6,7,8].  The  finite  element  method  is  especially  appropriate  for 
the  problem  of  penetration  of  electromagnetic  fields  into  an  inhomogeneous , 
absorbing  body.  To  treat  such  problems,  one  couples  the  finite  element  procedure 
inside  the  body  with  an  integral  equation  or  other  technique  on  the  surface  of  the 
absorbing  body.  For  some  work  on  this,  see  [9,10,11].  A  particularly  successful 
coupling  method  has  been  developed  by  Waterman  [12,13].  In  Waterman's  method,  an 
integral  equation  is  used  on  the  boundary,  and  the  external  fields  are  expanded  in 
a  series  of  harmonic  vector  fields.  This  eliminates  the  singularity  from  the 
kernel  of  the  equation  and  gives  a  representation  of  the  approximate  solution  in 
terms  of  fields  which  already  satisfy  the  differential  equations  of  the  problem. 

It  should  be  noted,  however,  that  although  the  resulting  Integrands  no  longer 
contain  singularities,  they  contain  oscillating  terms,  and  care  in  their  evaluation 
is  still  required.  Examples  of  further  work  along  this  line  are  contained  in 
[14,15,16]. 


5.  Goldstein,  C.,  "Numerical  Methods  for  Helmholtz  Type  Equations  in  Unbounded 
Domains,"  BNL-26543,  Brookhaven  Laboratory,  Brookhaven,  N.Y.,  1979. 

6.  McDonald,  B.  H.  and  Wexler,  A. ,  "Finite  Element  Solution  of  Unbounded  Field 
Problems,"  IEEE  Transactions  on  Microwave  Theory  and  Techniques,  Vol.  MTT-20, 
pp.  841-847,  1972. 

7.  Bettes,  P. ,  "Infinite  Elements,"  International  J.  for  Numerical  Methods  in 
Engineering,  Vol.  11,  pp.  53-64,  1977. 

8.  Kriegsmann,  G.  and  Morawetz,  C.  S. ,  "Numerical  Solutions  of  Exterior  Problems 
with  the  Reduced  Wave  Equation,"  Journal  of  Computational  Physics,  28, 

pp.  181-197,  1978. 

9.  Zienkiewicz,  0.  C.,  Kelly,  D.  W. ,  and  Bettess,  P.,  "The  Coupling  of  the  Finite 
Element  Method  and  Boundary  Solution  Procedures,"  International  J. 

for  Numer.  Methods  in  Engineering,  Vol.  11,  pp.  355-375,  1977. 

10.  Brezzl,  F.  and  Johnson,  C.,  '*0n  the  Coupling  of  Boundary  Integral  and  Finite 
Element  Methods,"  Report  77.15R,  Dept,  of  Computer  Science,  Chalmers 
University  of  Technology,  Gbteborg,  Sweden,  1977. 

11.  Brezzi,  F. ,  Johnson,  C.,  and  Nedelec,  J.  C. ,  "On  the  Coupling  of  Boundary 
Integral  and  Finite  Element  Methods,"  Report  39,  Centre  de  Mathematiques 
Appliquees,  Ecole  Poly technique,  Paris,  1978. 

12.  Waterman,  P.  C.,  "Matrix  Formulation  of  Electromagnetic  Scattering," 

Proc.  IEEE.  Vol.  53,  pp.  805-812,  1965. 

13.  Waterman,  P.  C. ,  "New  Formulation  of  Acoustic  Scattering,"  J.  of  the  Acoustical 
Soc. ,  Vol.  45,  pp.  1417-1429,  1969. 

14.  Barber,  P.  and  Yeh,  C.,  "Scattering  of  Electromagnetic  Waves  by  Arbitrarily 
Shaped  Dielectric  Bodies,"  Applied  Optics,  Vol.  14,  pp.  2864-2872,  1975. 

15.  Barber,  P.,  "Resonance  Electromagnetic  Absorption  by  Nonspherical  Dielectric 
Objects,"  IEEE  Transactions  on  Microwave  Theory  and  Techniques,  Vol.  MTT-25, 
pp.  373-381,  1977. 

16.  Morgan,  M.  A.  and  Mel,  K.  K. ,  "Finite  Element  Computation  of  Scattering  by 
Inhomogeneous  Penetrable  Bodies  of  Revolution,"  IEEE  Transactions  on  Antennas 
and  Propagation.  Vol.  AP-27,  pp.  202-214,  1979. 
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The  method  proposed  here  has  similarities  with  the  method  of  Waterman,  in 
that  we  also  use  harmonic  vector  fields.  Our  starting  point  is  different, 
however,  in  that  we  view  the  problem  as  a  system  of  differential  equations  coupled 
with  a  complicated  non-local  boundary  condition.  (See  [17]  for  a  similar  point  of 
view.)  The  method  was  first  developed  for  the  penetration  problem  [18]  and  was 
described  in  a  preliminary  way  in  [19]  and  [20]. 

2.  MATHEMATICAL  FORMULATION 


Let  0  be  a  bounded,  perfectly  conducting,  three  dimensional  domain 
with  boundary  r.  We  allow  the  possibility  that  0  is  a  disconnected  domain,  in 
which  case  r  is  a  collection  of  disjoint  closed  surfaces.  Let  Qg  be  the  set  of 
points  not  in  Q.  We  are  given  an  incident  electromagnetic  wave,  E  ,  H®,  of 
frequency  u.  Thus,  eP  and  H®  are  vector  fields  defined  In  all  space,  which 
satisfy  the  reduced  Maxwell's  equations 

7  *  E°  -  ctH°,  7  x  H°  -  BE0. 

If  the  wave  is  propagating  in  a  vacuum,  a  *  iwpg,  8  -  -iweg,  where  eg  - 
8.8574  x  10“12  meters/farad,  pg  -  1.2566  x  10“®  meters/henry.  In  general, 
a  -  lag,  B  ■  -iBg,  where  oig,  Bg  are  positive  real  numbers.  Let  E,  H  be  the 
electromagnetic  wave  resulting  from  the  scattering  of  the  incident  wave  by  the 
perfectly  conducting  body  fi.  Let  »  E  -  Ep,  H*  *  H  -  be  the  scattered  wave. 
Then  E(x),  H(x)  are  defined  for  xeQg,  and  are  determined  by  the  following  set  of 


equations: 

(2.1a) 

7  x  E  -  oH 

xeQo 

(2.1b) 

7.  x  H  ■  BE 

xen0 

(2.2) 

E  x  n  -  0, 

xer 

(2.3a) 

E1,  H1  -  0(r_1) , 

(2.3b) 

ir  »  S1  -  V*b  -1  ' 

o(r_1). 

r-x» 

(2.3c) 

^  x  h1  +  ^80/a0  E1  - 

o(r-1) , 

r-«» 

17.  MacCamy,  R.  C. ,  "Variational 

Procedures  for  a  Class  of  Exterior  Interface 

Problems,"  Report  79-9,  Department  of  Mathematics,  Carnegie-Mellon  University 
Pittsburgh,  Pa.,  1979. 

18.  Aziz,  A.  K.  and  Kellogg,  R.  B.,  "Finite  Element  Analysis  of  a 

Scattering  Problem,"  Report  BN-917,  Inst,  for  Phys.  Sci.  and  Tech.,  Univ.  of 
Md.,  College  Park,  Md.,  1979. 

19.  Aziz,  A.  K.  and  Kellogg,  R.  B.,  "A  Scattering  Problem  for  the  Helmholtz 
Equation,"  in  Advances  in  Computer  Methods  for  Partial  Differential  Equations- 
III.  R.  Vichnevetsky,  editor,  IMACS,  1979. 

20.  Kellogg,  R.  B.,  "A  Scattering  Problem  for  Maxwell's  Equations,"  Ibid. 
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Here  n  denotes  the  outward  pointing  unit  normal  to  r,  r  ■  |x|,  and  denotes 
the  unit  vector  in  the  radial  direction. 

The  system  (2.1a,b)  comprises  a  system  of  six  partial  differential  equations 
in  the  six  unknowns  consisting  of  the  components  of  £  and  H.  The  equations  (2.1), 
(2.2),  (2.3)  form  an  exterior  boundary  value  problem  for  this  system;  (2.2)  is  the 
boundary  condition  on  T,  and  (2.3)  are  the  boundary  conditions  at  infinity.  This 
problem  has  been  treated,  for  example,  in  [21] .  where  it  is  shown  that,  for 
reasonable  surfaces  T  and  incident  waves  E®,  H®,  the  problem  has  a  unique  solution. 

It  is  convenient  to  discuss  the  boundary  value  problem  in  terms  of  the 
scattered  wave.  For  this,  let  f  be  a  tangential  vector  field  on  T,  and  consider 
the  boundary  value  problem  defined  by 

(2.4a)  V  x  E1  -  oH1,  xeflg, 

(2.4b)  V  x  h1  -  gE1,  xc&q, 

(2.5)  n  x  *  f ,  xeT, 

and  (2.3).  From  [21],  it  is  known  that  this  problem  has  a  unique  solution.  If 
£  *  -n  x  E°,  then  the  solution  E1,  H1  of  (2.4),  (2.5),  (2.3)  gives  the  scattered 
wave  for  the  original  problem.  Thus,  it  suffices  to  solve  (2.4),  (2.5),  (2.3). 

Let  Jf  be  a  tangential  vector  field  on  r.  We  define  another  tangential  field, 
Kf ,  as  follows.  Let  E^-,  H*  be  the  solution  of  (2.4),  (2.5),  (2.3),  and  let 
Kf  *  n  «  H*.  The  operator  K  maps  tangential  vector  fields  into  tangential  vector 
fields.  We  also  define  a  bilinear  form,  B,  on  tangential  vector  fields.  If  f^  and 
£  are  two  tangential  vector  fields  on  T,  we  define 

B(f,  S)  m  |  "*I  x  K£dr  • 

r 

The  operator  K  and  the  bilinear  form  B  are  used  in  the  formulation  of  our 
numerical  method.  We  prove  that  the  bilinear  form  is  symmetric. 

Lemma  1.  B(f,  &)  •  B(g,  f) . 

Proof.  Let  be  the  solution  of  (2.4),  (2.5),  (2.3),  and  let  E^,  be 

the  solution  of  (2.4),  (2.5),  (2.3)  with  replaced  by  £.  Let  Br  be  a  ball  of 
radius  r  with  boundary  Sr.  Let  r  be  chosen  so  large  that  &CBr,  and  let  JIq  r  * 

&0  OBr.  We  have 

BE1  •  E2  -  oH1  •  H2  -  V^H1  x  E2)  . 


MUller,  C. ,  Foundations  of  the  Mathematical  Theory  of  Electromagnetic  Waves, 


Springer-Verlag,  New  York,  1969. 
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Integrating  this  over  ,  we  obtain 

u,  r 

n-H1  x  E2dr  +  j  e^H1  x  E2dr 

r  s 

r 

The  first  term  on  the  left  side  is 


[BE1*!2  -  aH1 »H2]dx. 


0,r 


|  n*E2  x  ^dr  -  j  n*£  x  Kfdr 


r 

-  b(*,  f) . 

Using  (2.3),  we  find  that  the  second  term  on  the  left  side  is 

j  I2*^.  x  J^dr-  -^Bq/Oq  j  E^E^dr  +  5(r) , 


where  6(r)  +  0  as  r  +  «,  Since  B(g.,  f)  is  independent  of  r,  it  follows  that  the 
expression 


[  (BE1*!2  -  aH1  *H2)dx  +  o  E^^dT 


0,r 


has  a  finite  limit  as  r  +  »,  and  that 

(2.6)  B(£,  f)  -  lim  J  (BE1*!!2  -  aH1*H2)dx  +  /Bq/oq  |  EX*E2dr 

*~'%T  S- 


11  2  2 

Since  the  right  side  of  (2.6)  is  unchanged  if  E.  ,  H  and  IJ  ,  H  are  interchanged, 
the  left  side  of  (2.6)  is  unchanged  if  f_  and  £  are  interchanged,  and  the  lemma  is 
proved. 


We  now  show  that  the  bilinear  form  is  definite.  This  will  be  used  in  the 
next  section  to  show  that  the  finite  system  of  equations  which  our  method 
produces  always  has  a  solution.  To  state  the  result,  we  let-z  denote  the  complex 
conjugate  of  a  complex  number  z. 


Lemma  2.  If  B(f_,  f)  ■  0,  then  _f  -  0. 


Proof.  Set  4  ■ 
(2.7)  ReB(f,  7)  - 


£_  in  (2.6),  and  note  that  oi  and  B  are  imaginary. 


(yy  iin 

r-*» 


Then 


In  particular,  we  find  that  the  limit  on  the  right  side  of  (2.7)  exists.  Suppose 
B(£,  f)  *  0.  Then 
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lim  d>  IE1 1 2dr  -  0. 


It  follows  from  a  theorem  of  Relllch  [21,  Theorem  15]  that  E1  =  0.  Hence  f.  *  0. 

3.  THE  APPROXIMATION  SCHEME 

We  describe  our  numerical  method  In  terms  of  the  scattered  wave  E^.  Let  £  be 
any  tangential  vector  field  on  r.  Then  since  n  x  Jjl  ■  -n  x  E®  on  r,  we  have 

n'E1  x  K£  -  n  x  E1*^ 

*  -n  x  E°-K& 

-  -n*E°  x  Kg^. 


Integrating  this  over  r,  and  letting  f_  denote  the  tangential  projection  of  the 
restriction  of  E*  to  r,  we  obtain 


B(f,  &) 


-♦ 


x  Kgdr. 


r 


This  equation  is  the  basis  for  our  numerical  method.  We  pick  a  finite  dimensional 
collection  S  of  tangential  vector  fields  on  T,  and  we  define  our  approximate 
solution  _f  eS  by 

(3.1)  B(f,  g)  -  -j  n-E°  x  K^dr,  gcS. 

r 


The  system  (3.1)  gives  rise  to  a  finite  system  of  linear  equations  whose  solution 
determines  the  vector  field  ftS.  This  approximation  scheme  seems  to  suffer  from 
two  defects.  It  is  not  clear  how  to  obtain  the  approximate  scattered  field  in  0q 
from  _f,  and  it  is  not  clear  how  to  evaluate  the  operator  K  which  appears  in  (3.1). 
These  defects  are  overcome  by  a  judicious  choice  of  subspace,  which  we  now  describe. 


Let  x*eft  be  given,  and  let  (r,0,$)  be  a  system  of  spherical  coordinates  with 
the  origin  at  x*.  In  Stratton  [22,  p.  416],  there  are  constructed  a  family  of 
vector  fields, 

(3.2)  ,1^  ,  m  -  0,1, •••,n,  n  -  1,2,*“. 

inn  mn 
o  o 

which  satisfy  the  following  properties. 

(1)  The  fields  are  regular  for  x  i  x*,  and  hence  are  regular  in 
(ii)  The  fields  satisfy 


7  x  n 


9 


22.  Stratton,  J.  A.,  Electromagnetic  Theory,  McGraw-Hill,  New  York,  1941. 
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1  x  Se  ’  • 

mn  mn 

o  o 

(iii)  The  fields  satisfy  (2.3). 

Note  that  we  must  take  the  Hankel  functions  zn(p)  *  h^(P)  In  Stratton’s  formulas 
to  satisfy  (iii). 

We  fix  an  integer  N  >  0,  and  let  Cfj  denote  the  collection  of  vector  fields 

(3.2)  for  0<m<n,  l<n<N.  Let  denote  the  collection  of  tangential  vector 

fields  on  T  which  are  of  the  form  n  *  (u|p),  where  UeCjj.  There  are  2N(N+2)  linearly 

independent  fields  in  Sjj.  If  &£%>  using  (ii),  we  may  easily  calculate  Kg.  If 

fe%  is  the  solution  of  (3.1),  then  f_  comes  from  a  vector  field  U  in  CN;  the  field 

11  is  the  desired  approximate  scattered  field,  and  may  be  easily  evaluated  at  points 
of  ftg.  We  have  therefore  shown  how  to  overcome  the  defects  of  using  (3.1). 

We  now  discuss  the  system  of  equations  arising  from  the  use  of  (3.1).  We 
arrange  the  fields  (3.2)  of  Cjj  in  a  definite  order  and  denote  them  by  Fp, 

1  <  V  i  M  •  2N(N+2).  We  let  «  n  *  (Ty|r).  Thus,  the  f^,  1  <  u  <  M,  are 

tangential  vector  fields  on  r  which  form  a  basis  for  Sty.  Writing  the  desired 
solution  f^  of  (3.1)  as  f  a  Ic^f^,  we  obtain  the  linear  system 

M 

(3.3)  2^.  tj  -  -  |  n  •  E°  x  Kf^dT,  1  i  u  5  M. 

v-i  r 

We  set  a  ■  B(f  ,  f  ),  and  we  let  A  ■  fa  ]  denote  the  coefficient  matrix.  We 
jjv  — y  —V  UV 

have  the  following 

Theorem.  The  complex  matrix  A  is  symmetric,  nonsingular,  has  nonsingular 
principal  minors,  and  admits  a  Cholesky  factorization  A  =  LLT. 

Proof.  The  symmetry  of  A  follows  from  Lemma  1.  The  nonsingularity  of  A,  and 
of  the  principal  minors  of  A,  follows  from  Lemma  2.  The  existence  of  the 
Cholesky  decomposition  then  follows  from  the  arguments  of  [23,  Thm.  3.1].  (Note, 
however,  that  A  is  complex  and  symmetric,  not  Hermit ian.) 

Remark  1.  The  right  side  of  (3.1)  could  be  written  in  terms  of  the  bilinear  form. 
However,  an  application  of  Lemma  1  would  not  enable  us  to  express  the  right  side 
in  terms  of  the  Incident  magnetic  field.  This  is  because  KEy  #  HO,  since  the 
incident  field  does  not  satisfy  (2.3).  To  avoid  this  possible  confusion,  we  have 
chosen  to  express  the  right  side  of  (3.1)  as  we  have  done. 

Remark  2.  An  error  analysis  for  the  method  proposed  here  will  appear  in  a  forth¬ 
coming  paper. 


23.  Stewart,  G.  W. ,  Introduction  to  Matrix  Computations,  Academic  Press,  New  York, 
1973. 
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4 .  COMPUTER  IMPLEMENTATION 


Two  computer  programs,  PCISH  (perfect  Conduction  In  Spherical  Harmonics)  and 
PSYM  (an  axisymmetric  version  of  PCISH),  have  been  developed  to  implement  the 
numerical  scheme  described  in  §3.  PCISH  is  the  more  general  of  the  two  since  it 
is  capable  of  computing  the  electric  (or  magnetic)  field  resulting  from  the 
scattering  of  an  arbitrary  electromagnetic  wave  E®  by  a  simply-connected, 
infinitely  conducting  body  0  of  arbitrary  shape.  PSYM  is  an  off-shoot  of  PCISH 
which  is  designed  to  handle  the  specific  class  of  problems  in  which  the  boundary 
of  ft  is  a  surface  of  revolution  and  the  incident  wave  is  a  plane  wave  propagat¬ 
ing  along  the  axis  of  symmetry. 


In  PCISH,  it  is  assumed  that  the  boundary  of  ft  can  be  suitably  approximated 
by  a  closed  surface  which  is  the  union  of  a  number  of  quadrilaterals.  The  vertices 
of  the  quadrilaterals  and  information  giving  the  assignments  of  vertices  to 
quadrilaterals  comprise  a  major  portion  of  the  input.  Other  input  parameters  are 
the  frequency,  orientation,  and  shape  of  the  incident  wave,  the  center  and  maximum 
order  N  of  the  subspace  C^,  the  quadrature  order,  and  the  coordinates  of  the 
points  in  space  at  which  the  scattered  field  is  to  be  calculated.  The  input  to 
PSYM  is  similar  (except  for  specifying  the  incident  wave)  but  is  much  simpler 
since  the  user  need  only  supply  a  suitable  number  of  points  lying  on  a  curve  T 
which  generates  the  boundary  of  ft  by  rotation  about  an  axis  of  symmetry.  The 
actual  rotation  of  T  is  carried  out  implicitly  by  the  program. 

Most  of  the  computations  performed  by  PCISH  and  PSYM  center  around  the  evalua¬ 
tion  of  the  coefficient  matrix  and  the  right-hand  side  of  the  linear  system  (3.3). 
In  view  of  the  Theorem,  both  PCISH  and  PSYM  use  Gauss  elimination  without  pivoting 
to  solve  the  matrix  equation. 

As  for  the  evaluation  of  the  coefficient  matrix  and  right-hand  side,  PCISH 
performs  a  two  dimensional  quadrature  (based  on  the  trapezoidal  rule)  over  each 
quadrilateral  of  the  surface  representing  the  boundary  of  ft.  PSYM,  on  the  other 
hand,  performs  only  a  one-dimensional  quadrature  (also  based  on  the  trapezoidal 
rule)  over  a  piecewise  linear  arc  approximating  T,  which  is  formed  by  joining  the 
given  input  points  on  T  by  line  segments.  The  reason  that  only  a  one  dimensional 
quadrature  is  necessary  is  that  in  the  case  in  which  3ft  is  a  surface  of  revolution, 
the  quadrature  in  the  $  direction  has  been  carried  out  exactly  by  hand  (requiring 
only  the  integration  of  some  trigonometric  polynomials)  and  the  corresponding 
formulas  placed  in  the  program.  In  both  PCISH  and  PSYM,  the  quadrature  order  is 
an  input  parameter. 


The  assumption  that  the  boundary  of  ft  is  a  surface  of  revolution  has  another 
large  advantage  in  that,  due  to  certain  symmetries  which  exist  in  this  case,  many 
of  the  matrix  entries  are  zero.  Moreover,  if  the  incident  wave  E®  is  assumed  to 
be  a  plane  wave  propagating  along  the  axis  of  symmetry,  certain  of  the  right-hand 
side  entries  also  become  zero,  and  the  originally  2N(N+2)-dimensional  system 
reduces  to  a  2N-dimensional  system  where  N  is  again  the  maximum  order  of  the 
subspace  Cjj.  Instead  of  using  all  of  the  vector  fields  (3.2),  we  may  take 
subspaces  generated  only  by  the  vector  fields 


m  . 

— oln 


Seln, 


The  advantage  of  such  a  reduction  is  obvious,  especially  for  large  values  of  N. 
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In  order  eo  evaluate  the  Integrand  In  the  bilinear  form  at  each  quadrature 
point,  one  needs  to  be  able  to  calculate  the  Hankel  functions  h^  and  the 
associated  Legendre  polynomials  P®  used  in  the  definition  of  the  fields  of  (3.2). 
The  routines  used  by  PCISH  and  PSYM  are  as  follows. 

The  Hankel  function?'  ■  jn  +  iyn  are  computed  at  a  given  point  p  >  0  using 
a  backward  recursion  formula  to  calculate  jn,  the  nth  order  spherical  Bessel 
function  of  the  first  kind,  and  a  forward  recursion  formula  to  calculate  yn,  the 
nth  order  spherical  Bessel  function  of  the  second  kind.  Starting  with  Xjj_^  *  1 
and  x„  <■  0  for  M  sufficiently  large,  we  recursively  compute  x„,  n  *  0,1, •••,M-2 
by 

/_  .  xn+l 

xn  -  (2n  +  3)  —  -  xn+2  . 


We  then  set 


j  (p)  ,  Stop.  x 

nvv/  pxQ  n 

As  for  the  yn,  we  set  yQ(p)  *  -  -^cosp,  y^(p)  *  -  ^cosp  -  ^sinp  and  compute 


yn<»> 


2n-l 


Vl(p)  •  yn-2(p) 


We  then  set  l\|(p)  ■  jn<P)  +  iyn(p)»  One  also  needs  certain  derivatives  of  the 
h^  and  for  these,  we  employ  the  following  useful  formula  (see  10.1.21  of  [24]) 


P  Q 


h1  .  -  -  h1 
n-1  p  n 


To  calculate  the  associated  Legendre  polynomials  Pn  of  degree  n,  order  m,  at 
a  given  point  n»-  1  -  n  -  1  we  first  observe  that  using  the  Rodriguez  formula  for 

pn  “  pn»  we  have  „ 

.  def  2  f  A.W 

P  (n)  -  (1  -  n  )  - - — 

n 


dn 


2.2  f  1 

"  "  dn”  L2”"1 


-  U 

„  (1  -  n*)*  d1 


,  2  .  .n 

(n  -  l) 


dn 


o 

2.2  ,m+n 


2nnl 


dn 


m+n 


,  2  ,.n 

(n  -  1) 


24.  Abramowitz,  M.  and  Stegun,  I.  A. ,  Handbook  of  Mathematical  Functions. 

U.S.  Government  Printing  Office,  1964. 
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Therefore, 


m 

,,  2.2  ,2m  „ 

p>>  -  t1- ^  <*2  -  d" 

TO  .TO  .  ,  ZTO 

2  m!  dn 


(1  -  n2)2(2m) 1 
2®m! 


TO-1 


(2»H2»rn  I2(S-1?IL-  (1  .  n2)  2  >^77 

2m  2®”1(m-l) ! 


-  (2m-l)  Vl  -  n2  pj£j(n>  . 

Since  Pg(n)  =  1,  we  may  recursively  calculate  P®(n)  for  any  m.  To  obtain  P®(n) 
for  m  4  n,  we  observe  that  Pg(n)  =  0  for  n  <  m  and  apply  the  recursion  formula 
(see  8.5.3  of  [24]) 


(n  -  «H)P^fl(n)  “  (2n+l)nP®(n)  -  (tt**)P®_1(n)  • 

We  also  need  the  derivatives  ^  P®.  Using  the  fact  that  P®(n)  =  0  for  all  m, 
we  apply  the  formula  (see  8.5.4  of  [24]) 

dp® 

(1  -  n2)  dT  "  (n4lB)Pn-l  "  nTlPn 

for  n  t  1,  m  *  0. 


Once  the  linear  system  has  been  solved  by  Gauss  elimination,  the  approximate 
scattered  field  JE-*-  is  assembled  and  may  be  evaluated  at  any  given  point  x  in 
space.  The  relevant  quantity  for  many  applications  is  the  radar  cross-section 
(RCS)  defined  by 


(4.1) 


o(x) 


ll°(x)|2 


where  E  is  the  incident  wave  and  R  is  the  distance  of  the  point  x  from  the  body 
Q.  To  calculate  the  scattered  field  at  infinity  (i.e.  the  farfield)  one  takes 
the  limit  of  (4.1)  asR->®.  Computationally,  this  is  achieved  by  using  an  asymptotic 
form  of  the  Hankel  functions  when  the  scattered  field  E^is  assembled  from  the 
solution  of  the  linear  system. 


5.  NUMERICAL  RESULTS 


Both  PCISH  and  PSYM  have  been  developed  and  executed  on  the  CDC  6500 
computer  located  at  NSWC/White  Oak.  These  programs  have  been  applied  to  a 
number  of  relatively  simple  model  problems  in  order  to  test  the  program 


14 


NSWC  TR  80-245 


capabilities  and  to  compare  computed  results  with  results  found  in  the  literature. 
Since  all  of  the  problems  attempted  to  date  have  been  axially  symmetric,  the  more 
specialized  program  PSYM  was  used  to  obtain  the  numerical  results  presented  below. 

As  mentioned  in  14,  the  maximum  order  N  and  center  x*  of  the  subspace  Cjj  and 
the  quadrature  order  are  input  parameters.  In  order  to  verify  our  results  for  a 
given  body  ft,  we  manipulate  these  quantities  in  the  following  way.  We  first 
determine  the  proper  value  for  N  by  successively  increasing  this  parameter  until 
the  RCS  stabilizes  to  some  pre-determined  number  of  significant  digits.  However, 
since  an  increase  in  N  causes  an  increase  in  the  oscillation  of  the  associated 
Legendre  polynomials  used  in  their  definition,  we  must  simultaneously  increase 
the  quadrature  order  to  insure  an  accurate  evaluation  of  the  matrix  system  (3.3). 
For  a  given  value  of  N,  the  quadrature  order  is  increased  until  the  RCS  again 
stabilizes  to  a  pre-determined  tolerance.  Once  these  quantities  have  been 
determined,  we  perform  a  consistency  check  by  varying  the  subspace  center  x*. 

Since  this  should  in  theory  have  no  effect  on  the  farfleld  profile,  this  seems  to 
be  a  fairly  rigorous  test  of  the  computed  values.  Obviously,  for  each  new  body 
ft  this  process  involves  a  fair  number  of  program  runs,  but  it  is  hoped  that  a 
further  investigation  of  our  method  will  lead  to  a  better  understanding  of  the 
Interplay  among  these  various  quantities  and  their  relationship  to  the  other 
relevant  parameters  of  the  problem. 

In  the  case  in  which  ft  is  an  infinitely  conducting  sphere,  an  exact  solution 
is  known  and  can  be  found  in  many  classical  texts  on  electromagnetic  theory 
(see  e.g.  [22]).  PSYM  produced  extremely  good  farfield  backscattering  results 
in  this  case,  even  well  into  the  resonance  region,  as  can  be  seen  by  comparing 
Fig.  1  with  a  similar  plot  obtained  from  the  exact  solution  on  p.  148  of  [25]. 

The  case  in  which  ft  is  a  right  circular  cone  is  also  studied  in  the  litera¬ 
ture  [25].  For  our  tests,  we  used  a  right  circular  cone  circumscribed  about  a 
1  m.  sphere  with  a  15®  half-angle  at  the  vertex.  The  plane  wave  E®  was  incident 
upon  the  vertex  and  propagated  along  the  axis  of  symmetry.  Fig.  2  shows  the 
farfleld  backscattering  results  produced  by  PSYM  in  the  Rayleigh  region  and 
extending  a  short  way  into  the  resonance  region.  These  results  seem  to  agree 
with  other  computed  and  experimental  results  reported  in  [25]  (in  particular, 
see  the  figure  on  p.  392  of  [25]). 

For  some  applications,  it  is  desired  to  calculate  the  scattered  field  close 
to  the  body  ft.  This  presents  no  difficulties  for  PSYM,  and  Fig.  3  displays  some 
partial  results  in  this  direction  for  the  same  conical  scatterer  as  described 
above.  For  each  of  the  frequencies  1MHz,  10MHz,  and  50MHz  (corresponding  to 

•  .027,  .273,  and  1.365  respectively  in  Fig.  2),  the  scattered  field  was 

evaluated  on  circles  of  radii  4m. ,  10m. ,  and  50m.  lying  in  the  plane  parallel  to 
E°  which  contains  the  axis  of  symmetry  of  the  cone.  In  order  to  obtain  a  farfield 
profile,  calculations  were  also  made  on  a  fourth  circle,  concentric  with  the  other 


25.  Ruck,  G.  T. ,  editor.  Radar  Cross  Section  Handbook.  Vol.  1,  Plenum  Press, 
New  York,  1970. 
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three  circles,  which  had  an  effectively  infinite  radius.  The  frequency  1  MHz  is 
in  the  Rayleigh  region  and  yields  a  constant  RCS  in  the  farfield.  The  frequency 
50  MHz  is  in  the  resonance  region,  and  the  frequency  10  MHz  is  on  the  border 
between  the  Rayleigh  and  resonance  regions  (see  (26]). 


26. 


Crispin,  J.  W. ,  Jr.,  and  Maffett,  A.  L. ,  "Radar  Cross-Section  Estimation  for 
Simple  Shapes,"  Proc.  IEEE.  Vol.  53,  833-847,  1965. 
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